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ABSTRACT 



Moments of inertia were experimentally determined and longitudinal and 
lateral/directional static and dynamic stability and control derivatives were estimated for a 
fixed wing Unmanned Air Vehicle (UAV). Dynamic responses to various inputs were 
predicted based upon the estimated derivatives. A divergent spiral mode was revealed, 
but no particularly hazardous dynamics were predicted. The aircraft was then 
instrumented with an airspeed indicator, which when combined with the ability to 
determine elevator deflection through trim setting on the flight control transmitter, 
allowed for the determination of the aircraft's neutral point through flight test. The neutral 
point determined experimentally corresponded well to the theoretical neutral point. 
However, further flight testing with improved instrumentation is planned to raise the 
confidence level in the neutral point location. Further flight testing will also include 
dynamic studies in order to refine the estimated stability and control derivatives. 
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I. INTRODUCTION 



A. MISSION NEED 

Unmanned Aerial Vehicles (UAVs) are gaining acceptance as an integral part of the 
operations of today's armed forces. Preceding and during Operation Desert Storm, UAVs 
flew a variety of missions including reconnaissance, targeting for gunfire support, and 
battle damage assessment. Desert Storm provided the first-ever combat test of Unmanned 
Aerial Vehicles by U.S. forces [Ref. 1]. Reviews of lessons learned from Desert Shield 
and Desert Storm reveal the outstanding successes of UAVs. 

There are many examples of effective UAV use during the war. In one instance, a 
commander of a Marine task force was able to monitor UAV imagery of Kuwait as his 
task force approached the city, revealing the exact reaction of the Iraqi forces to Marine 
armor, artillery, and troop movements. The Navy used UAVs to search for mines, spot 
for gunfire support, perform reconnaissance missions for SEAL teams, and search for 
Iraqi Silkworm sites, command and control bunkers, and anti-aircraft artillery sites. The 
Marines quickly reacted to an Iraqi attack into Saudi Arabia observed by a UAV and 
decisively crushed the Iraqi invasion with airborne Cobras and Harriers. The Army 
provided their Apache pilots with route reconnaissance acquired from UAVs shortly 
before the Apache missions. Spotting for air strikes and naval gunfire support 



1 



became so successful that Iraqi soldiers were seen attempting to surrender to UAVs as 
they flew overhead [Ref. 2], 

UAVs have many advantages over manned aircraft which help account for their 
effectiveness. First, the cost of a UAV is a very small fraction of the cost of a manned 
aircraft. The Pioneer, for example, costs approximately $500,000 for the aircraft and 
$500,000 for the onboard camera, bringing the total cost of the package to a mere one to 
two percent of the cost of most manned tactical aircraft. UAVs are also extremely flexible 
with respect to launch platform. They can be launched from small fields, truck beds, and 
practically any ship in the Navy inventory. UAVs are frequently very hard to detect with 
radar or infrared systems due to their small size, composite construction, small engines 
(sometimes electric motors), and their slow speed. UAVs also have an advantage from 
being unmanned. The aircraft is not limited by the "g" tolerance of a pilot or by pilot 
fatigue. Finally, the best advantage of all is that when a UAV crashes or is lost to enemy 
fire, there is no search and rescue mission required, no prisoners of war taken, and no loss 
of life. 

UAVs are not without their problems, however. Acquisition and support programs 
are relatively new and underdeveloped. This results in a UAV force that is relatively small 
in number of aircraft, and small and inexperienced in terms of personnel. The Pioneer 
showed signs of these underlying problems during Desert Storm. Six Pioneers were 
damaged badly enough to require return to the factory for repair due to no intermediate 
level maintenance facilities being available. Five Pioneers were lost due to mechanical 
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malfunction or operator error [Ref. 2], Due to the small number of available aircraft, 
losses are very costly and operator errors need to be avoided if possible. Thorough and 
frequent training through simulation can keep operators performing at their peak. 

B. STATEMENT OF OBJECTIVE 

To provide realistic simulation for training, accurate aerodynamic data for the 
aircraft must be available. Aerodynamic parameters of the aircraft can be determined from 
its physical characteristics, wind tunnel tests, and flight tests of actual or scale models. It 
is the objective of this work to provide the ground work for determining whether flight 
test can effectively be used to accurately estimate the static and dynamic stability and 
control characteristics of a fixed wing UAV. Flying-qualities parameters were estimated 
using analytic techniques, and the resultant flight dynamics were simulated to provide 
expected behavior for future test flights. 
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n. BACKGROUND 



A. EVOLUTION OF AN AIRCRAFT 

For any aircraft to get from an idea to an actual flying vehicle, properties such as 
stability and control derivatives, which determine the flying characteristics of the aircraft, 
must be determined. These derivatives are used to size control surfaces, design flight 
control systems, and program training devices such as simulators. There are typically 
three ways to determine or estimate derivatives. 

The first method of determining an aircraft's derivatives begins early and continues 
throughout the design process. This step involves determining derivatives mathematically 
from physical characteristics of the aircraft. Requiring little more than a few 
well-educated engineers and some calculators or desktop computers, this method is 
relatively simple. Derivatives can be reasonably approximated, but must be refined 
through other methods. 

The second method of determining derivatives is through aerodynamic 
measurements of wind tunnel testing. This method is more complicated than simple pencil 
and paper calculations due to the necessity for model making and wind tunnel operation, 
but much better results can be achieved. Derivatives must still be refined, however, due to 
factors such as scale effects and interference from wind tunnel walls and 
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supporting hardware. Also, dynamic effects are often difficult to properly account for in 
most wind tunnels. 

The final step approach to determining the derivatives of an aircraft is through flight 
test. This method is by far the most expensive and complicated way to determine the 
aircraft's derivatives, but it is also the most precise and complete. Scaled flight test 
provides a viable option for this third method. 

B. PIONEER UAV 

In June 1982, Israeli forces very successfully used UAVs as a key element in their 
attack on Syria. Scout and Mastiff UAVs were used to locate and classify SAM and AAA 
weaponry and to act as decoys for other aircraft. This action resulted in heavy Syrian 
losses and minimal Israeli losses. A year and a half later, the U S. Navy launched strikes 
against Syrian forces in the same area with losses much higher for the Navy than those of 
the Israelis [Ref. 3], 

The Commandant of the Marine Corps, General P. X. Kelly, recognized the 
effectiveness of the Israeli UAVs. Secretary of the Navy John Lehman then initiated 
development of a UAV program for the U.S. Anxious to get UAVs to the fleet. Secretary 
Lehman stipulated that UAV technology would be off-the-shelf [Ref. 4], After the 
contract award to AAI Corporation of Baltimore, Maryland for the Pioneer UAV, 
developmental and operational testing took place concurrently. This approach resulted in 
quick integration of the Pioneer into the fleet. Unfortunately, such quick integration into 
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the fleet can result in problems identified during operational use which had not been fully 
explored in the test and evaluation process. 

The UAV Office at the Pacific Missile Test Center (PMTC, now the Naval Air 
Warfare Center, NAWC, Weapons Division, Pt. Mugu) was tasked with Developmental 
Test and Evaluation of the Pioneer. Testing revealed the following concerns which 

warranted further investigation [Refs. 5,6]: 

1 . discrepancies in predicted with flight-tested rate of climb, time to climb, 
and fuel flow at altitude; 

2. apparent autopilot-related pitch instability; 

3. tail boom structural failure; 

4. severely limited lateral control; 

5. slow pitch response causing degraded maneuverability at high gross 
weights; 

6. insufficient testing to determine the effects of the new wing on flight 
endurance. 

The Target Simulation Laboratory at Pt. Mugu was tasked to develop a computer 
simulation of the Pioneer in order to provide cost-effective training for pilots. 

Aerodynamic data were needed to provide the stability and control derivatives necessary 
for the simulation as well as to answer questions concerning basic flying qualities of the 
Pioneer. 

In order to provide support to the research being done at Pt. Mugu and to provide 
for future UAV project support, a research program was begun at the Naval Postgraduate 
School (NPS). An instrumented half-scale radio-controlled model of the Pioneer was used 
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for the research at NPS. Research performed included wind tunnel tests, flight tests, and 
numerical modeling. 

Initial NPS research on the Pioneer, performed by Capt. Daniel Lyons, involved a 
computer analysis of the Pioneer in its original configuration and with a proposed larger 
tail. A low order panel method (PMARC) was used for the aerodynamic analysis. Static 
longitudinal and directional stability derivatives, the neutral point, and crosswind 

limitations were calculated. Drag polars were constructed using the component buildup 

/ 

method for profile drag, and drag reduction measures were considered [Ref. 7], 

In conjunction with Capt. Lyons work, Lt. James Tanner conducted wind tunnel 
tests to determine propeller efficiencies and thrust coefficients for drag studies [Ref. 8], 

Lt. Tanner also conducted flight tests to determine power required curves and drag polars 
[Ref. 8], Capt. Robert Bray later conducted wind tunnel tests of a 0.4-scale model at 
Wichita State University to determine static stability and control derivatives [Ref. 9], 
Aerodynamic data obtained by Capt. Lyons and Bray have been supplied to PMTC to be 
used for simulation. 

Lt. Jim Salmons performed initial flying qualities flight testing using an onboard data 
recording system in order to determine static stability parameters. Unfortunately, 
vibration problems with the onboard recorder rendered much of the data unusable [Ref. 
10 ]. 

Following up on Lt. Salmons' work, Lt. Kent Aitcheson installed the CHOW-1G 
telemetry system, designed by Lt. Kevin Wilhelm, in an attempt to alleviate the vibration 
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problem experienced by Lt. Salmons. The new flight test configuration was used to test 
static longitudinal and lateral-directional stability characteristics of the Pioneer. The 
vibration problem experienced by Lt. Salmons was overcome, though not enough data 
were acquired for a complete and thorough analysis of the Pioneer's characteristics. Much 
insight was gained, however, concerning instrumentation. Resolution needed to be 
improved for flight control position indication [Refs. 1 1,12]. 

Lt. Paul Koch conducted further flight tests of the Pioneer with the CHOW-1G 
telemetry system. Static longitudinal stability results from the flight tests correlated well 
with theoretical predictions and with simulations of a full-scale Pioneer. Electromagnetic 
interference with the flight control system at the test site resulted in loss of the half-scale 
model Pioneer before further data collection and analysis could be performed [Ref. 13]. 

C. PARAMETER ESTIMATION 

Parameter estimation is used to derive stability and control derivatives from dynamic 
flight test data. Lcdr. Robert Graham successfully used the Modified Maximum 
Likelihood Estimation (MMLE) technique with MATLAB software to analyze simulated 
and actual flight test data of several aircraft. Analysis of simulated UAV data revealed the 
effect of signal-to-noise ratio on the estimator, and the need for proper control-surface 
excitation for a particular response [Ref. 14], Cdr. Patrick Quinn analyzed flight test data 
from the Marine Corps BQM-147 UAV using both MMLE and a more robust non-linear 
model, pEst, and compared the results of the two approaches. Noise was found to be a 
problem when the system response was in the same general frequency as the noise. 
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Limited available data and noise at the frequency expected for the system's response 
prevented a successful resolution of all stability and control derivatives of interest. The 
pEst model was found to be the estimator of choice when aircraft maneuvers exceed what 
is generally considered reasonable for linear approximation of flight dynamics [Ref. 15], 
While previous work with UAVs at NPS has at times been frustrating, much has 
been learned, especially concerning the most challenging method of aircraft analysis, flight 
test. This work initiates the implementation of the lessons learned from previous work 
into dynamic simulation and flight testing of a generic UAV in order to properly prepare 
the UAV lab at NPS for further support of future projects and to demonstrate the value of 
scaled UAV flight test to the fleet. 



9 



m. THE AIRCRAFT 



A. GENERAL DESCRIPTION 

The "Bluebird", shown in Figures 1-3, is a high-wing tricycle-gear radio-controlled 
airplane. It is constructed of wood, foam, composites, and metal. It is powered by a 
Sachs-Dolmar 3.7-cubic-inch two-stroke gasoline engine which drives a 24 inch 
two-bladed wood propeller. It is controlled by a nine-channel pulse-code-modulated 
Futaba radio operating at 72.710 MHz. To enhance reliability, the Bluebird has two 
receivers which share control of the aircraft. The left receiver controls the left aileron, 
elevator, and flap, and engine ignition and onboard electronics package cut-off, while the 
right receiver controls the right aileron, elevator, and flap, and rudder, nose-wheel 
steering, and throttle. The Bluebird can fly within visual range for approximately 1.5 
hours. Table 1 describes physical specifications of the Bluebird. 

B. STABILITY DERIVATIVES 

The initial estimates of the stability derivatives of the Bluebird were made using the 
physical characteristics of the aircraft such as airfoil data , geometric measurements, 
relative positions of aircraft components, mass, and weight [Refs. 16-19]. The assumed 
flight condition with the associated aircraft configuration is described in Table 2. 
Nondimensional stability and control derivatives estimated are shown in Table 3, and 
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dimensional stability and control derivatives estimated are shown in Table 4. MATLAB 



programs written for the physical and derivative calculations appear in Appendix A. 

TABLE 1 

SPECIFICATIONS 



Length 

Height 

Wing Airfoil (est.) 

Horizontal Stab. Airfoil (est.) 

S„ing (SJ 

S, 

S v 



c . • 

b .. 

b, • 

by • 

AR 

AR, 

ARy 

V H 

V v - 

V' . 



.... 9.84ft. 
.... 3.04ft. 
... GO 769 
NACA 4412 
.. 22.38 ft. 2 
. . 4.701 ft. 2 
.. 1.277 ft. 2 
... 1.802 ft. 
... 1.281 ft. 
... 12.42 ft. 
.... 3.67ft. 
.... 1.21 ft. 

6.89 

2.86 

1.14 

.... 0.6274 
.... 0.0247 
.... 0.0023 
. . . 5.381 ft. 
. . . 5.381 ft. 



TABLE 2 

FLIGHT CONDITION/AIRCRAFT CONFIGURATION 



Weight 



ZZ 

Velocity 

Altitude 

Density 

Center of Gravity 
r 

^Ltrim 

A o A, rim 

Elevator (rim 



57.79 lbs. 

.... 12.58 slug*ft 2 
.... 13.21 slug*ft 2 

19.99 slug*ft 2 

. 60 mph/88 ft/sec 

800 ft MSL 

0.002327 slugs/ft 3 

27% of m.a.c. 

0.2866 

3.8 degrees 

1.6 degrees 
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Figure 1 
Side View 





Front View 
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TABLE 3 

NONDIMENSIONAL DERIVATIVES 



C L 


0.2866 


Cd 


0.0358 


C U 


4.1417 


Cdo 


0.0311 


Cl 


1.5787 


Cd a 


0.1370 


C„\ 


-1.0636 


Cm 


-4.6790 


Cl, 


3.9173 


Cm] 


-11.6918 


Clu 


0.4130 


Cd v 


0.0650 


Cmu 


-1.2242 


c, P 


-0.3100 


Cn$ 


0.0484 


C/ p 


-0.0330 


Cy P 


0.0000 


Cn p 


-0.0358 


Ci, 


-0.3579 


Cy, 


0.0967 


Cn 


-0.0526 


Cl 


0.0755 


Cysa 


0.0000 


Cnso 


-0.0258 


Clu 


0.2652 


Cy&r 


........ 0.0697 


Cns, 


-0.0326 


C/„ 


0.0028 


Cd g 


0.0000 


Cy, 


0.0000 



TABLE 4 

DIMENSIONAL DERIVATIVES 



X u -0.0914 /sec 

X&e -7.2961 ft/sec 2 

Z a -468.9852 ft/sec 2 

Z q 4.5027 ft/sec 

M u 0.0000 /ft*sec 

-1.3178 /sec 

Mse -33.6730 /sec 2 

Y p 0.0000 ft/sec 

Ysa 0.0000 ft/sec 2 

Lp -6.5787 /sec 2 

L r 1.0613 /sec 

Lsr 0.5589 /sec 2 

N p -0.3167 /sec 

Nsa -3.2375 /sec 2 



Z a 16.7894 ft/sec 

Z u -0.7312 /sec 

Z* 1.8146 ft/sec 

Z& -46.368 ft/sec 2 

M a -29.2559 /sec 2 

M q -3.2928 /sec 

Yp -34.8021 ft/sec 2 

Y r 0.7663 ft/sec 

Y&r -7.8282 ft/sec 2 

L p -5.0281 /sec 

Lsa 52.7966 /sec 2 

Np 6.0593 /sec 2 

N r -0.4647 /sec 

N S r -4.0900 /sec 2 
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C. MOMENTS OF INERTIA 



Having accurate moments of inertia is critical to ensuring the accurate prediction of 
aircraft dynamics. Direct calculation of a model's moments of inertia by consideration of 
the contributions made by individual parts is impractical and inaccurate. Determination of 
the moments of inertia by test is much more practical and precise. The changes in a 
model's moments of inertia due to addition or subtraction of equipment or structure can be 
calculated directly, thereafter. 

In the determination of moments of inertia by test, the aircraft is hung from the 
ceiling and swung. Using the period exhibited and the principles of compound pendulums, 
the moments of inertia of the model can be extracted [Ref. 20-21], In order to calculate 
the moment of inertia about all axes, the model must be hung from the ceiling and swung 
three different ways, each such that as the aircraft swings, it is rotating about the axis of 
interest. The Bluebird was hung by chain and swung as pictured in Figures 4-6. 
Specifications for the geometry of each test can be found in Appendix B. 

Reference 21 provides equation (1) for calculating the moment of inertia of a 
swinging model. 



t w s z s _ 2 WmZm 

1 - 4 „- Pm+s - rPs - g 



(i) 



W is the weight, Z is the distance from pivot to center of gravity, p is the period, and g is 
the gravitational constant. M and S are subscripts designating either the model or the 

support. 
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Figure 4 

1^ Test (not to scale) 
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I w Test (not to scale) 
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I a Test (not to scale) 
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It was determined that swinging the support (the chains) in the configuration it 
would be in when supporting the model would not be possible, since the chains would not 
maintain their positions without the model in place. Equation (1) was therefore 
manipulated in order to treat the chains as long slender rods and to calculate their 
moments of inertia as such [Ref. 22], 



The new form of equation (1) is 



j _ W „ 2 " 

1 ~ 4n J Pm+s s l 



w M zh 



3g 



( 2 ) 



where Ls is the length of a chain and the summation is taken over all chains (four in this 
case). In particular, there were two "long" chains and two "short" chains, all having a 
weight per unit length © . Appropriate substitutions were made to yield 



r _ + 2 o>(Lshort + Llong)]Zm+s n 2m f r 2 . 7 - 2 

1 ~ ^ 'M+S 1 3 SHORT + L LONG) V J 

Having the equation in this form fixes the values for all variables except 
Z M+ s, Z M , and Pm+s . These three variables were measured for each of the three 
configurations and calculations were made. Four periods were timed during the swing 
tests, and the tests were repeated ten times. The moments of inertia calculated are shown 
in Table 2. 
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IV. PREDICTED DYNAMICS 



A. GENERAL 

When flight testing an aircraft, it is important to attempt to predict the results of the 
testing prior to the actual flights. The information provided by the predictions can be used 
in briefing the pilot as to what to expect from the aircraft and also to avoid any potentially 
dangerous flight regimes. Longitudinal and lateral-directional dynamics of the Bluebird 
have been predicted using the stability and control derivatives estimated in chapter three 
along with computational methods based upon the six equations of motion as described in 
references 18 and 23. Computational programs were written in MATLAB and appear in 
Appendix C. 

B. LONGITUDINAL DYNAMICS 

The MATLAB program named "longnat.m" uses the full (4x4) longitudinal plant 
and the MATLAB "impulse" function to determine the short and long period natural 
responses. A diary of the values computed and output by the program provides 
eigenvalues, damping ratios, and damped and undamped natural frequencies for both 
modes as shown in Table 5. 
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TABLE 5 

LONGITUDINAL DATA 



Eigenvalues 
Damping ratio 

Undamped Natural Frequency 
Damped Natural Frequency 



Short Period Long Period 
-5.083 +/- 4.861 i -0.037 +/- 0.400i 



0.723 

7.03 rad/sec 
4.86 rad/sec 



0.093 

0.401 rad/sec 
0.399 rad/sec 



Figures 7 and 8 show the combined short and long period natural response. It can be seen 
that the short period response is almost completely damped out after only two seconds. 
Response beyond two seconds is primarily long period. The phugoid mode is seen to be 
lightly damped. 

The short and long period response to a unit step elevator input is determined by 
using the full (4x4) longitudinal plant and the MATLAB "step" function in the MATLAB 
program named "stepper. m". Figures 9 and 10 show the short and long period response 
to a step input. In the first plot, the long period can be seen to be much more heavily 
excited than the short period mode. Again, it can be seen that the short period response is 
almost completely damped out after about two seconds. Due to held elevator input, the 
long-term response is to trim to a new angle of attack, while pitch rate dies out. 

The MATLAB program named "n_step.m" uses the full (4x4) longitudinal plant and 
the MATLAB "step" function to determine the normal acceleration or load factor 
response to a unit step elevator input. Figure 1 1 shows the short period normal 
acceleration response to a unit step input. It should be noted that normal acceleration is 
defined opposite in sign to load factor. The long-period response is seen to be much 
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Figure 7 

Natural Response (short time) 
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Figure 8 

Natural Response (long time) 
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Figure 9 

Response to Step Input (short time) 




Figure 10 

Response to Step Input (long time) 
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more critical than the short-period response in generating large load factors for this 
aircraft. 




Figure 1 1 

Normal Acceleration Response to Step Input 



The longitudinal response to given initial conditions is determined using the full 
(4x4) longitudinal plant and the MATLAB "initial" function in the MATLAB program 
named "homogen.m". The initial conditions ( u/U = .34, alpha = 5 deg, q = 8.8 deg/sec, 
theta = -0.8 deg) are provided for all states from the step input results after initial 
dynamics have died out (approximately 15 seconds). Figures 12 and 13 show the 
longitudinal response to the initial conditions. The response in angle of attack is small 
while a significant pitch rate is developed in both the short-period and long-period 
responses. The long-period is again seen to be lightly damped. 

The MATLAB program named "doublet. m" uses the full (4x4) longitudinal plant 
and the transition matrix (using the MATLAB "exponential" function) to find the 
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Figure 12 

Homogeneous Response (short time) 
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Figure 13 

Homogeneous Response (long time) 
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response to a unit pitch doublet input at approximately the short period damped natural 
frequency. Figure 14 shows the input pitch doublet, and Figure 15 shows the response. 
Though the doublet is commonly used to observe only the short-period response, the 
angle of attack response indicates the long-period response is also excited. 

The response (transfer function gain and phase) to a harmonic elevator input is 
found using the full (4x4) longitudinal plant and the MATLAB "bode" function in the 
MATLAB program named "sp_bode.m". Figures 16 and 17 show the gains and phase 
shifts from an elevator frequency sweep. The angle of attack gain can be seen decreasing 
from a peak at a very low excitation frequency; the short-period response is masked by the 
long-period gain. In actuality, the low-frequency response is of little interest. Pitch rate 
shows a maximum gain at approximately 2.4 radians/second corresponding to an elevator 
cycle period of approximately 2.6 seconds. 

The MATLAB program named "n_bode.m" uses the full (4x4) longitudinal plant 
and the MATLAB "bode" function to find the normal acceleration or load factor response 
(transfer function gain and phase) to a harmonic input. Figures 18 and 19 show the 
normal acceleration gains and phase shifts from an elevator frequency sweep. The short 
and long-period damped natural frequencies can be observed at 4.86 radians/second and 
0.399 radians/second respectively where their respective gains peak. While the gain 
demonstrated at the long-period damped natural frequency is quite significant, the aircraft 
would not be expected to be excited at that frequency. An excitation near the 
short-period damped natural frequency is much more likely. Flight test pilots and 
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Figure 14 

Unit Elevator Doublet Input 




Figure 15 
Doublet Response 
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Frequency, rad/sec 



Figure 16 

Frequency Response (Gain) 




Figure 17 

Frequency Response (Phase) 
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Figure 18 

Normal Acceleration Response to a Harmonic Input (Gain) 




Normal Acceleration Response to a Harmonic Input (Phase) 
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engineers should be aware of the fact that excitation at that frequency will produce 
approximately 0. 15 "g's" per degree of elevator deflection or that it takes a harmonic 
amplitude of about 6.7 degrees of elevator to produce one "g" of acceleration. The 
Bluebird has a maximum of 15 degrees of up elevator and 12 degrees of down elevator 
available. This would equate to approximately -0.8 to +3.2 "g's", which is easily tolerated 
structurally. 

C. LATERAL-DIRECTIONAL DYNAMICS 

Undamped natural frequency, damped natural frequency, damped natural period, 
and damping ratio for the Dutch roll mode; time constant for the roll mode; and time 
constant and time to double or half for the spiral mode are calculated by the MATLAB 
program named "lat_dir.m" using the full (4x4) lateral-directional plant. Values calculated 
are described as follows. 



Dutch roll mode: 



Undamped natural frequency 


2.65 rad/sec 


Damped natural frequency 


2.62 rad/sec 


Damped natural period 


2.40 sec 


Damping ratio 


0.148 


Roll mode: 


Time constant 


0. 195 sec 


Spiral mode: 


Time constant 


-29.28 sec 


Time to double amplitude 


20.29 sec 



The Dutch roll mode is observed to be lightly damped. Also, the spiral mode is divergent 
with a fairly short time to double bank angle. 
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The MATLAB program named "rudkick.m" uses the full (4x4) lateral-directional 
plant and the MATLAB "impulse” function to find the response to a unit rudder impulse. 
The program also finds the bank angle at the end of 100 seconds which is approximately 
four degrees in the direction of the rudder kick. Figure 20 shows the unit rudder kick 
response in bank and sideslip angles. Sideslip lags bank angle by approximately 0.6 
seconds or 80 degrees of phase and is approximately 60 percent larger in magnitude. 

The response to an initial sideslip is determined by the MATLAB program named 
"sideslip.m" using the full (4x4) lateral-directional plant and the MATLAB "initial" 
function. Initial conditions for sideslip and bank angle are provided by a steady-sideslip 
condition where bank angle is ten degrees, and sideslip angle is 13.7 degrees. Figure 21 
shows the homogeneous response to the initial sideslip. The lightly damped Dutch roll 
mode can be seen by observing the oscillations of both bank angle and sideslip angle, while 
the divergent spiral mode can be observed by the increasing bank angle. 

The MATLAB program named "roll.m" uses the full (4x4) lateral-directional plant 
and the MATLAB "bode" function to find the roll angle response (transfer function gain 
and phase) due to a harmonic aileron input. Figures 22 and 23 show the roll angle gains 
and phase shifts from the harmonic aileron input. The maximum high-frequency roll angle 
gain is seen to occur at approximately 2.8 radians/second which corresponds to an aileron 
input period of approximately 2.2 seconds. 
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Figure 20 

Unit Rudder Kick Response 




Figure 21 

Homogeneous Response to an Initial Sideslip 
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Frequency, rad/sec 



Figure 22 

Roll Angle Gain due to a Harmonic Aileron Input 




Figure 23 

Roll Angle Phase Shift due to a Harmonic Aileron Input 
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D. REMARKS 



Analysis has been conducted for the most common modes. The Bluebird 
demonstrates no particularly dangerous flying qualities. Particular characteristics of the 
Bluebird's behavior worth noting include a very heavily damped short period longitudinal 
mode, a lightly damped phugoid mode, a lightly damped Dutch roll mode, and a divergent 
spiral mode. Additional analysis could be easily done by analogy to those modes all ready 
analyzed. Further analysis might include, for example, the response to a step aileron 
input, aileron doublet, or elevator impulse ("stick rap"). 
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V. FLIGHTIEST 



A. TEST DESCRIPTION 

The relative positions of an aircraft's center of gravity and neutral point are critical 
to the aircraft's longitudinal stability and handling qualities. In order for a conventional 
(aft tail) aircraft to be statically stable, its center of gravity must be forward of its neutral 
point. The farther the center of gravity is in front of the neutral point (relative to the 
chord length), the more statically stable the aircraft is. As the center of gravity is moved 
aft beyond the neutral point, the aircraft becomes statically unstable and more 
maneuverable. 

As an aircraft's center of gravity is moved aft toward the neutral point, less and less 
change in elevator trim is required to achieve steady, level flight for a given change in 
airspeed. This fact can be used to experimentally determine the aircraft's neutral point. 
Reference 23 provides equation (4) which describes the relationship between change in 
pitching moment coefficient with change in coefficient of lift and the required change in 
elevator deflection with change in coefficient of lift. 

(4) 

As the center of gravity moves aft and approaches the neutral point, dC m /dC L approaches 
zero. Since V H and C LSe are constants, d5e/dC L must also approach zero. 
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In order to determine the neutral point by flight test, the aircraft must be flown at 
various centers of gravity. At each particular center of gravity, the aircraft is trimmed at 
several airspeeds (coefficients of lift). Each airspeed and the corresponding elevator 
deflection are recorded. Plots are then generated showing elevator deflection versus 
coefficient of lift at each center of gravity tested. A line is then best fit through the data 
points for each center of gravity. The slope of this line represents d5e/dC L for that 
particular center of gravity. Finally, d5e/dC L versus center of gravity is plotted. A best fit 
line is then drawn through these data points. Using the best fit line just drawn, the center 
of gravity where d5e/dC L equals zero is the aircraft’s neutral point. 

B. INSTRUMENTATION 

In order to acquire the airspeed of the Bluebird, a simple, commercially-available 
airspeed indicator was installed. The Digicon TT-01 Tele Tachometer/ASI senses the 
spinning of a small wind-driven propeller blade using a cadmium disulphide optical sensor. 
The frequency of the changes in light intensity sensed is transmitted to a hand-held 
receiver which converts the frequency to airspeed which can be read directly in real-time. 
The manufacturer's claimed accuracy is +/- 0.5 feet per second. The airspeed indicator 
was installed on the left wing of the Bluebird by mounting it on the end of a boom which 
extended approximately 18 inches (approximately 80 percent of c-bar) in front of the 
leading edge of the wing. 

Measuring elevator trim was done in an indirect manner. The flight control 
transmitter has a digital display which can show elevator trim on a generic scale of -100 to 
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100. Prior to flight test, elevator deflection in degrees was measured for various trim 
settings, and the relationship between transmitter displayed trim number and actual 
degrees of elevator deflection was determined. This calibration allowed for the 
determination of elevator deflection in degrees during post flight analysis. 

C. TEST PROCEDURES 

In order to fly the Bluebird at various centers of gravity, an access panel in the aft 
fuselage of the aircraft was modified to accept added weights. The aircraft was then 
configured with various amounts of weight and the location of the center of gravity 
determined for each configuration. Centers of gravity were determined by weighing each 
wheel with a strain-gage balance. 

A total of seven flights were flown at various centers of gravity. Flights were kept 
short in order to minimize the shift in the center of gravity due to fuel bum. The location 
of the center of gravity was determined both before and after the flight, and the average 
center of gravity was used for calculations. Shifts in center of gravity during testing were 
kept to one percent or less of the wing chord. 

Each flight consisted of 12 passes at various airspeeds. Airspeed and trim setting 
for each pass were recorded for post flight analysis. Additionally, air pressure, air 
temperature, and aircraft weight were noted in order to calculate coefficient of lift. 

Post flight analysis of the data was conducted in accordance with the procedures 
described in Section A, Test Description. The resulting graphs are shown in Figures 
24-31. 
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Degrees of Elevator 
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25 FEB 94, Flight 1 
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Figure 24 
Flight 1 of 7 
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Figure 25 
Flight 2 of 7 
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C.G. = .3198 

04 MAR 94, Flight 2 
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Figure 27 
Flight 4 of 7 
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C.G. = .3493 

04 MAR 94, Flight 3 
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Figure 28 
Flight 5 of 7 
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Figure 29 
Flight 6 of 7 
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Degrees of Elevator 


C.G. — .3848 

04 MAR 94, Flight 5 
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Figure 30 
Flight 7 of 7 
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D. RESULTS 



The neutral point estimated by flight test was found to be about 54 percent of the 
mean aerodynamic chord. The neutral point estimated through conventional calculations 
was approximately 52.8 percent of the mean aerodynamic chord. 

While such a close match between the neutral point locations found by each method 
is highly encouraging, one must not overlook the scatter in the data presented in Figure 
3 1 . Assuming a normal distribution of the data, there is a 68 percent confidence that the 
experimentally determined neutral point is within five percent of the wing mean 
aerodynamic chord of the location determined experimentally. Such data scatter is 
unacceptable for accuracies required, and the tests will be repeated when the improved 
instrumentation under development comes on line. 

For the flight testing done, two sources of potential error are of particular interest. 
First, it was very difficult to ensure a perfectly level pass of the aircraft at each particular 
airspeed. The pilot was positioned on the ground and the plane was flying approximately 
one hundred feet above him. This is not an optimum vantage point for the pilot. Ideally, 
the pilot would like to be elevated (such as in a tower) such that the aircraft is at eye level 
for each pass. Future plans include using an altitude-hold autopilot to maintain the desired 
trim conditions. Also, due to wind gusts and possibly small unintentional changes in 
aircraft attitude, the airspeed was not always steady, and therefore somewhat difficult to 
read consistently. This second problem will be overcome by the improved data acquisition 
system. 
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VL SUMMARY 



A determination of the moments of inertia by the component contribution method 
was decided to be too cumbersome for analysis of the complete aircraft. Moments of 
inertia were found through a compound pendulum analysis, and appear reasonable. Future 
changes to the configuration of the aircraft can be accounted for by considering the 
contribution of the component added, removed, or relocated. 

Initial estimates of stability and control derivatives were made by conventional 
aircraft-design-type methods. Such methods involve considering characteristics of the 
aircraft such as lift-curve slopes of airfoils and the locations of particular aircraft 
components relative to each other. Programs were written in MATLAB in such a manner 
that future changes to the aircraft or different flight conditions can be accounted for 
quickly and accurately. The initial estimates of the derivatives can be used in the 
preliminary design of a flight controller. 

The initial estimates of the stability and control derivatives were used to predict the 
dynamic response of the aircraft to several different excitations. MATLAB programs 
were written to predict the dynamics, and any future modifications to the aircraft or flight 
conditions can be accounted for easily. Several characteristics of the aircraft's responses 
are worth noting. First, the short period longitudinal response was found to be heavily 
damped. Little or no evidence of the short period response can be observed beyond 
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approximately two seconds. Second, the long period longitudinal response was found to 
be lightly damped. The long period response can still be observed relatively easily after 
100 seconds. Similarly, the Dutch roll mode is also fairly lightly damped and can be 
observed easily for ten to 15 seconds. Finally, the spiral mode was found to be divergent 
with a time to double bank angle of approximately 20 seconds. While this mode is not 
particularly dangerous for a radio controlled aircraft which is flown visually (open loop), 
the pilot should be aware of this tendency of the aircraft in order to avoid trim problems. 
It is also recommended that flight controller designs incorporate bank angle feedback for 
wing leveling. 

The neutral point of the aircraft was found to be 53 to 54 percent of the wing mean 
aerodynamic chord by two different methods. The data collected through flight test were 
somewhat scattered; it is estimated that the neutral point location is known within +/- five 
percent of the wing mean aerodynamic chord with a confidence of about 68 percent. 
Further flight testing with improved instrumentation is recommended to raise the 
confidence of the neutral point estimation. It is recommended that flight tests with 
improved instrumentation continue to verify or update all predicted stability and control 
derivatives. 
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APPENDIX A: MATLAB PROGRAMS USED TO ESTIMATE 
STABILITY AND CONTROL DERIVATIVES 

1. blbrdfcl.m 

% Eric J. Watkiss 
% AA0810 Thesis 

% File for Bluebird data which change with flight condition 
% blbrdfcl.m 

% Last Update: 02 MAR 94 



g = 32.174; 
Wlmg = 24.02 
Wrmg = 23.91 
Wng = 11.07 


% Acceleration due to gravity 
% Weight on left main in lbs 
% Weight on right main in lbs 
% Weight on nose gear in lbs 


Umph = 60 
Ufps = 88.0 


% Flight speed in miles per hour 
% Flight speed in feet per second 


rho = .002327 


% Air density in slugs/(cubic ft) 


Ixx = 12.58 
Iyy = 13.21 
Izz = 19.99 
Ixz = 0 


% Moment of inertia about x-axis 
% Moment of inertia about y-axis 
% Moment of inertia about z-axis 
% Assumed! ! !••!'!!!!!!!!!!!!!!!! 


LD = 8; 


% Lift to drag ratio 


thetanaut = 0; 


% Initial pitch angle 
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2. bluebird. m 



% Eric J. Watkiss 
% AA0810 Thesis 

% File for Bluebird data which are fixed 
% bluebird.m 

% Last Update: 02 MAR 94 



ac = .479; 


% 


ai = 3; 
% 


% 


alphaOl = -6. 5* pi/ 180; 
% 

% 


% 


b = 12.42; 


% 


bt = 3.67; 


% 


bv = 1.208; 


% 


cbar = 1.802; 
% 


% 


CLalphaafw = 5.443; 
% 

% 


% 


CLalphaaft = 5.587; 
% 

% 


% 


CLalphaafv = 2*pi; 
% 

% 


% 


CMac = -.06; 
% 


% 


ct = 1.281; 


% 


c4tail = 7.878; 
% 

% 


% 


c4wing = 2.497; 
% 


% 


daOdde = .625; 
% 

% 


% 


daOddr = .675; 
% 

% 


% 


deda = .4; 


% 



Aileron chord in ft. 
distance from centerline to 
inner edge of aileron in ft. 
a.o.a. for zero lift (radians)ao = 6; 
distance from centerline to 
outer edge of aileron in ft. 

Span of wing in ft 
Span of horizontal tail in ft. 

Height of vertical tail in ft. 

Mean aerodynamic chord (m.a.c.) 
in feet 

Lift curve slope of wing 
airfoil (GO 769) in per 
radian 

Lift curve slope of horizontal 
tail airfoil (NACA 4412) in per 
radian 

Lift curve slope of vertical 
tail airfoil (flat plate) in per 
radian 

Coefficient of moment about 
aero. ctr. (GO 769) 
m.a.c. of horizontal tail in ft. 
Location of quarter chord of 
horizontal tail in feet from 
firewall 

Location of quarter chord of 
wing in feet from firewall 
Section flap effectiveness 
for 33% flap (elevator) 

Abbott and Doenhoff p. 190 
Section flap effectiveness 
for 38% flap (rudder) 

Abbott and Doenhoff p. 190 
Downwash angle derivative 
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% 




Df = 1; 


% 


eO = 0; 


% 


ee = .8; 


% 


g = 32.174; 


% 


hac = .245; 


% 


% 




it = 4.83 *pi/l 80; 


% 


lewing = 2.047; 


% 


% 




letail = 7.557; 


% 


% 




% 




mg = 37.595/12; 


% 


% 




ng = . 75/12; 


% 


% 




S = 22.380; 


% 


Sr = .547; 


% 


St = 4.701; 


% 


Sv = 1.277; 


% 


Wf = .67; 


% 


ybar = b/4; 


% 


zv = .5; 


% 


% 




Zwf = .5; 


% 



% 



estimated from Perkins/Hage 
Depth of fuse, in ft. 

Assumed epsilon naught 
Assumed span efficiency factor 
gravitational constant 
Location in percent chord of 
aero. ctr. (NACA 4412) 
Incidence angle of hor. tail 
Location of leading edge of wing 
in feet from firewall 
Location of leading edge of 
horizontal tail in feet from 
firewall 

Location of main gear in ft 
from firewall 

Location of nose gear in ft 
from firewall 

Reference (wing) area in sq. ft. 
Rudder area in sq. ft. 

Horizontal tail area in sq. ft. 
Vertical tail area in sq. ft. 

Width of fuse, in ft. 

Spanwise location of m.a.c. 

Vert, tail height to m.a.c. 
(estimated) 

Verticle height of wing 
above fuse. C.L. in ft. 
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3. dderiv.m 



% Eric J. Watkiss 
% AA0810 Thesis 

% File to calculate dimensional derivatives 
% dderiv.m 

% Last Update: 12 FEB 94 

% Run nondimensional derivative program 
ndderiv 

% Calculate dynamic pressure 

qbar = ,5*rho*Ufps A 2; % ft lbs 

Malpha = CMalpha*qbar*S*cbar/Iyy; % per second A 2 

Mq = CMq*(cbar/(2*Ufps))*qbar*S*cbar/Iyy; 

% per second 

Malphadot = CMalphadot*(cbar/(2*Ufps))*qbar*S*cbar/Iyy; 

% per second 

Xu = -2*CD*qbar*S/(m*Ufps); % per second 

Zu = -2*CL*qbar*S/(m*Ufps); % per second 

Zalphadot = CLalphadot*(cbar/(2*Ufps))*(qbar*S/m); 

% ft per second 

% ft per second 

% per ft second 

% ft per second A 2 

% ft per second A 2 

% per second A 2 

% ft per second A 2 

% ft per second A 2 



Zq = CLq*(cbar/(2*Ufps))*(qbar*S/m); 
Mu = 0; 

Xde = -l*CDde*qbar*S/m; 

Zde = -l*CLdelta*qbar*S/m; 

Mde = CMde*qbar*S*cbar/Iyy; 

Xalpha = (CL - CDalpha)*qbar*S/m; 
Zalpha = -l*(CLalphaw+CD)*qbar*S/m; 
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YB = CyB*qbar*S/m; 


% ft/second A 2 


LB = ClB*qbar*S*b/Ixx; 


% l/second A 2 


NB = CnB*qbar*S*b/Izz; 


% 1 /second A 2 


Yp = Cyp*b*qbar*S/(2*Ufps*m); 


% ft/sec 


Yr = Cyr*b*qbar*S/(2*Ufps*m); 


% ft/sec 


Lp = Clp*(b/(2*Ufps))*qbar*S*b/Ixx; 


% 1/sec 


Np = Cnp*(b/(2*Ufps))*qbar*S*b/Izz; 


% 1/sec 


Lr = Clr*(b/(2*Ufps))*qbar*S*b/Ixx; 


% 1/sec 


Nr = Cnr*(b/(2*Ufps))*qbar*S*b/Izz; 


% 1/sec 


Ydr = -l*Cydr*qbar*S/m; 


% ft/sec A 2 


Yda = 0; 


% ft/sec A 2 


Ldr = Cldr*qbar*S*b/Ixx; 


% l/sec A 2 


Lda = Clda*qbar*S*b/Ixx; 


% l/sec A 2 


Ndr = Cndr*qbar*S*b/Izz; 


% l/sec A 2 


Nda = Cnda*qbar*S*b/Izz; 


% l/sec A 2 



Malphaprime = Malpha + Malphadot*(Zalpha/Ufps); 
Mqprime = Mq + Malphadot; 

Mdeprime = Mde + Malphadot* (Zde/Ufps); 
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4. godallas.m 

% Eric J. Watkiss 
% AA0810 Thesis 

% File to get values of dimensional and nondimensional derivatives 
% godallas.m 

% Last Update: 04 FEB 94 

!del bluebird. dia 
diary bluebird. dia 
diary on 

% Run programs to calculate derivatives 
dderiv 

% Nondimensional Derivatives 

CL 

CD 

CDO 

CLalphaw 

CMalpha 

CDalpha 

CLalphadot 

CMalphadot 

CLq 

CMq 

CLdelta 

CDde 

CMde 

CyB 

CnB 

C1B 

Cyp 

Cnp 

Clp 

Cyr 

Cnr 

Clr 

Cyda 

Cnda 

Clda 

Cydr 

Cndr 



49 



Cldr 

CDq 

Cyq 

% Longitudinal Dimensional Derivatives 
Xu 

Xalpha 

Xde 

Zu 

Zalpha 

Zalphadot 

Zq 

Zde 

Mu 

Malpha 

Malphadot 

Mq 

Mde 

% Lateral/Directional Dimensional Derivatives 

YB 

Yp 

Yr 

Yda 

Ydr 

LB 

Lp 

Lr 

Lda 

Ldr 

NB 

Np 

Nr 

Nda 

Ndr 

diary off 
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5. ndderiv.m 



% Eric J. Watkiss 
% AA0810 Thesis 

% File to calculate nondimensional derivatives 
% ndderiv.m 

% Last Update: 12 FEB 94 

% Load Bluebird data with flight condition 
physcalc 

% Calculate coefficients of lift and drag 
CL = W/(.5*rho*Ufps A 2*S); 

CD = CL/LD; 

% Calculate lift curve slope of wing in per radian 
CLalphaw = CLalphaafw/(l+CLalphaafw/(pi*ee*AR)); 

% Calculate lift curve slope of horizontal tail in per radian 
CLalphat = CLalphaaft/(l+CLalphaaft/(pi*ee*ARt)); 

% Calculate lift curve slope of vertical tail in per radian 
CLalphav = CLalphaafv/(l+CLalphaafV/(pi*ee*ARv)); 

% Calculate change in hor. tail lift with change in elevator 
dcLtdde = daOdde * CLalphat; % per radian 

% Calculate change in vert, tail lift with change in rudder 
dcLvddr = daOddr * CLalphav; % per radian 

% Calculate zero lift pitching moment 
CMO = CMac + VH * CLalphat * (it + eO); 

% Calculate CMalpha in per radian 

CMalpha = CLalphaw*((h-hac)-VH*(CLalphat/CLalphaw)*(l-deda)); 

% Calculate change in aircraft lift with change in elevator 
CLdelta = dcLtdde*(St/S); % per radian 

% Calculate chng in aircraft pitching moment w. chng in elevator 
CMde = -1 *VH*dcLtdde; % per radian 

% Calculate angle of attack and elevator angle for trimmed flight 
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% 

% CM = CMO + CMalpha*alpha + CMde*de 
% Cl = CLalphaw*alpha + CLdelta*de 

% 

% 

% | CLalphaw CLdelta | | alpha | | CL | 

% I II 1=1 I 

% |_CMalpha CMde | |_dej |_-CMOJ 

% 

% A * X = C 

% 

A = [ CLalphaw CLdelta 
CMalpha CMde ]; 

C = [ CL 
-l*CMO ]; 

X = inv(A)*C; 

atrim = X( 1,1); % trim a.o.a. in radians 

etrim = X(2, 1), % trim elevator in radians 

% Calculate change in yawing moment with change in rudder 
% "rudder power" 

% assumes VF/Vinfinity = 1 

Cndr = - 1 * W*dcLvddr; % in per radian 

% Calculate CnB contribution from vert, tail 
% CnB = CLalphav*W*(W/Vinfinity) A 2*(l-dsigma/dbeta) 

% assumes VF/Vinfinity = 1 and dsigma/dbeta = 0 
CnB = CLalphav* W;% in per radian 

% Calculate change in rolling moment with change in sideslip 

% First calculate dihedral contribution from wing 
% Raymer p. 439 

CIBwf = -1.2* sqrt( AR) *Zwf‘ (Df+ W f)/b A 2 ; 

CIBw = ClBwCL*CL+ClBwf; 

% Next calculate contribution from fin 

% CIBv = -l*Clalphav*Vprime*(VF/Vinfinity) A 2*(l-dsigma/dbeta) 
% Assume VF/Vinfinity = 1 and dsigma/dbeta = 0 
CIBv = -l*CLalphav*Vprime; % in per rad 

% Combine CIBg and CIBv into C1B 
C1B = CIBw + CIBv; % in per rad 
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% Calculate "aileron power", Clda 
% See Smetana pp. 139-141 
Cldatau = Cldatauo - Cldataui; 

Clda = Cldatau*tau; % in per radian 

% Calculate change in yawing moment w. aileron deflection 
Cnda = 2*K*CL*Clda; % in per radian 

% Calculate side force due to yaw 

% By Smetana p. 107 

CyB = -.31; % in per radian 

% Calculate side force due to rudder 
Cydr = CLalphav*taur*Sv/S; % in per radian 

% Calculate side force due to aileron 
% By Smetana, p. 138 
Cyda = 0; 

% Calculate rolling moment due to rudder 
Cldr = Cydr*zv/b; % in per radian 

% Calculate change in drag due to change in elevator 

% Smetana pp. 95-100 

% Using Figure 26 at 8 degrees aoa 

CDde = ((. 1 55-.047)/(20*pi/l 80))*St/S; % in per radian 

% Calculate change in drag with change in aoa 

% Smetana pp. 64-65 

% Assuming dCDO/dalpha is negligible 

CDalpha = 2*CL*CLalphaw/(pi*ee*AR); % in per radian 

% Calculate change in pitching moment w.r.t. alphadot 
% Smetana pp. 78-81, etat assumed = 1 
CMalphadot = -2*CLalphat*deda*(ltprime/cbar)*... 
(lt/cbar)*(St/S); % in per radian 

% Calculate change in lift with pitch rate 
% Smetana pp. 82-85 

% Neglecting wing contribution, assuming etat = 1 
CLq = 2*(lt/cbar)*CLalphat*(St/S); % in per radian 
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% Calculate change in lift with alphadot 
% Smetana pp. 75-76 

CLalphadot = - 1 *CMalphadot/(lt/cbar); % in per radian 

% Calculate change in pitching moment w. pitch rate 
% Smetana pp. 87-88 
% Assuming etat = 1 

CMq = -2*(cbar/4-h)*abs(cbar/4-h)*CLalphaw/(cbar A 2) - ... 
2*(lt/cbar) A 2*CLalphat*(St/S); % in per radian 

% Calculate roll damping 

% Smetana pp. 122-125 

% Neglecting contribution from vertical tail 

Clp = -,475*(AR+4)/(2*pi*AR/CLalphaw+4); % in per radian 

% Calculate change in yawing moment due to rolling 
% Smetana pp. 126-129 
% Neglecting contribution from vertical tail 
Cnp = -1 *CL/8; % in per radian 

% Calculate change in side force with yaw rate 
% From Schmidt p. 3-23 
% Assume etat = 1 

Cyr = 2* W*CLalphav; % in per radian 

% Calculate change in rolling moment w. yaw rate 

% Schmidt p. 3-24 

% Tail contribution 

Clrv = (zv/b)*Cyr; % in per radian 

% Wing contribution 

Clrw = CL/4; % in per radian 

% Total 

Clr = Clrv + Clrw; % in per radian 

% Calculate yaw damping 

% Schmidt p. 3-25 

% Tail contribution 

Cnrv = -l*(lv/b)*Cyr; % in per radian 

% Wing contribution from Smetana p. 136 

CDO = CD-CL A 2/(pi*ee*AR); 

Cnrw = -.02*CL A 2-.3*CD0; % in per radian 
% Total 

Cnr = Cnrv + Cnrw; % in per radian 
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% The following 3 derivatives are negligible and taken to be 0 
CDq = 0; % in per radian 

Cyq = 0; % in per radian 

Cyp = 0; % in per radian 

% A few misc. calculations 

% Static Margin/Neutral Point 
statmar = CMalpha/(-l *CLalphaw) 
hn = statmar + h 
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6. physcalc.m 

% Eric Watkiss 
% AA0810 Thesis 

% File to calculate physical considerations 
% physcalc.m 

% Last Update: 04 FEB 94 

% Load fixed Bluebird data 
bluebird 

% Load flight condition 
blbrdfc 1 

% Calculate aircraft weight 
W = Wlmg + Wrmg + Wng; 

% Calculate aircraft mass 
m = W/g; 



% Calculate aspect ratio of wing 
AR = b A 2/S; 

% Calculate aspect ratio of hor. tail 
ARt = bt A 2/St; 

% Calculate aspect ratio of vert, tail 
ARv = bv A 2/Sv; 

% Calculate longitudinal center of gravity 
h = ((ng*Wng + mg * (Wlmg+W rmg))AV -lewing)/cbar ; 

% Calculate "tail length" from e g. to horizontal tail a.c. 
% same for horizontal and vertical 
It = c4tail - (lewing + h*cbar); 
lv = It; 

% Calculate "tail length" from c/4 wing to c/4 tail 
ltprime = c4tail - c4wing; 

% Calculate hor. tail volume coefficient 
VH = lt*St/(S*cbar); 
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% Calculate vert, tail volume coefficient (yaw) 

W = lv*Sv/(b*S); 

% Calculate vert, tail volume coefficient (roll) 

Vprime = zv*Sv/(b*S); 

% Unit antisymmetrical angle of attack for outer and inner 
% edge of aileron (See Smetana p. 141) 
antisymo = ao/(b/2); 

Cldatauo = .74; 
antisymi = ai/(b/2); 

Cldataui = .23; 
cacw = ac/cbar; 
tau - .52; 

% for yawing moment due to aileron, see p. 142, Smetana 
eta = ai/(b/2); 

K = -.17; 

% for side force due to rudder deflection, see Smetana p. 145 
vratio = Sr/Sv; 
taur = .62; 

% for rolling moment due to sideslip, See Raymer, Fig. 16.21, p. 439 
CIBwCL = -.04; 
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APPENDIX B: MOMENT OF INERTIA CALCULATIONS 



For the formula 

j + 2<o(L short + L long)]? m+s n M’.uZl/ 2<of r 2 , r 2 

1 ~ ^2 r Ar+S g 3 J SHORT + L LONG 



the following data were used for extraction of the moments of inertia. 



Fixed values: 



Weight of the model, 

Weight per unit length of chain. 

Length of short chain. 

Length of long chain. 

Gravitational constant, 

(adjusted for latitude and elevation) 


W M = 58.45 lbs 
co = 0.06133 lbs/ft 
L SHORT — 13.00 ft 
Llong — 15.00 ft 
g = 32.1472 ft/sec : 


Variable values: 






Distance from pivot to center of gravity 
of model and support. 

Distance from pivot to center of gravity 
of model. 

Average period of model and support. 


Zm+s = 12.95 ft 

Zm — 13.33 ft 
Pm+s = 4.109 sec 




Distance from pivot to center of gravity 
of model and support. 

Distance from pivot to center of gravity 
of model. 

Average period of model and support, 


Z M+S = 12.01 ft 

Z M = 12.35 ft 
Pm+s = 3.976 sec 




Distance from pivot to center of gravity 
of model and support. 

Distance from pivot to center of gravity 
of model, 

Average period of model and support. 


Z M + S = 12.49 ft 

Z M = 12.81 ft 
Pm+s = 4.077 sec 
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APPENDIX C: MATLAB PROGRAMS USED TO ESTIMATE 
AIRCRAFT DYNAMICS 



A. LONGITUDINAL DYNAMICS 
1. longnat.m 

% Determines the short and long period natural response 
% using the full longitudinal plant 

clear 

dderiv 

!del longnat.dia 
diary longnat.dia 

% Inertial Matrix 

in = [Ufps 0 0 0; 

0 Ufps-Zalphadot 0 0; 

0 -l*Malphadot 1 0; 

0 0 0 1 ]; 

% Aircraft Matrix 

an = [Ufps*Xu Xalpha 0 - 1 *g*cos(thetanaut); 

Ufps*Zu Zalpha Ufps+Zq -l*g*sin(thetanaut); 
Ufps*Mu Malpha Mq 0; 

0 0 1 0 ]; 

a = inv(in)*an; 
b = [1; 1; 1; 1]; 
c = eye(size(a)); 
d = [0; 0; 0; 0]; 

t = 0:01:100; 

[y,x,t] = impulse(a,b,c,d,l,t); 
clg 

figure(l) 

%plotting alpha and q, short period 
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plot(t,y(:,2)),axis([0 15-1.5 1 5]);grid;gtext('alpha’) 
hold on 

plot(t,y(:,3));gtext(’q’) 

xlabel(’Time, sec');ylabel(’Response Gain') 

%titleCNatural Response') 

hold off 

pause 

figure(2) 

%plotting u/U and theta, phugoid 
plot(t,y(:,l));axis([0 100 -2 2]);grid;gtext(’u/U') 
hold on 

plot(t,y( : ,4));gtext(’theta’) 
xlabel('Time, sec');ylabel('Response Gain') 
%titleCNatural Response'); 
hold off 

damp(eig(a)) 

diary off 
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2. stepper.m 



% Determines the short and long period step response 
% using the full longitudinal plant 

clear 

dderiv 

% Inertial Matrix 

in = [Ufps 0 0 0; 

0 Ufps-Zalphadot 0 0; 

0 -l*Malphadot 1 0; 

0 0 0 1 ]; 

% Aircraft Matrix 

an = [Ufps* Xu Xalpha 0 -l*g*cos(thetanaut); 
Ufps*Zu Zalpha Ufps+Zq -l*g*sin(thetanaut); 
Ufj)s*Mu Malpha Mq 0; 

0 0 1 0 ]; 

% Control Matrix 

bn = [Xde Zde Mde 0]'; 

a = inv(in)*an; 
b = inv(in)*bn; 
c = eye(size(a)); 
d =[0 0 0 0]'; 

t = 0:. 01:100; 

[y,x,t] = step(a,b,c,d,l,t); 

clg 

figure(l) 

%plotting alpha and q, short period 

plot(t,y(:,2)) 

axis([0 15 -5 3]) 

grid;gtext('alpha') 

hold on 

plot(t,y(:,3));gtext('q') 

xlabel('Time, sec');ylabelCResponse Gain') 

%title('Short Period Step Response') 
hold off 
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pause 



figure(2) 

%plotting u/U and theta, phugoid 

plot(t,y(:,l)) 

axis 

grid;gtext('u/U') 
hold on 

plot(t,y(:,4));gtext('theta') 
xlabel('Time, sec');ylabel('Response Gain') 
%title('Long Period Step Response'); 
hold off 
pause 

figure(3) 

%plotting alpha and q, long period 
plot(t,y(:,2)) 
grid;gtext('alpha') 
hold on 

plot(t,y(:,3));gtext('q') 
xlabel('Time, sec');ylabel('Response Gain') 
%title('Long Period Step Response') 
hold off 
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3. n_step.m 



% This program uses STEP function to determine 
% normal acceleration response to a unit step input 

clear 

dderiv 

% Inertial Matrix 

in = [Ufps 0 0 0; 

0 Ufps-Zalphadot 0 0; 

0 -l*Malphadot 1 0; 

0 0 0 1 ]; 

% Aircraft Matrix 

an = [Ufps* Xu Xalpha 0 -l*g*cos(thetanaut); 
Ufps*Zu Zalpha Ufps+Zq -l*g*sin(thetanaut); 
Ufps*Mu Malpha Mq 0; 

0 0 1 0 ]; 

% Control Matrix 

bn = [Xde Zde Mde 0]'; 

a = inv(in)*an; 
b = inv(in)*bn; 
c = [ Zu Zalpha Zq 0 ]; 
d = Zde; 

t =0:0.05:15; 

[y,x,t] = step(a,b,c,d,l,t); 

for i = 1 :length(t) 
y2 = y/32.174*pi/180; 
end 

% plotting g's/deg 
plot(t,y2);grid 

xlabel('Time, sec');ylabelCNormal Acceleration: g s/deg') 
%title(' Acceleration Response to Unit Step') 
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4. homogen. m 



% Determines the short and long period homogeneous response 
% using the full longitudinal plant 

clear 

dderiv 

% Inertial Matrix 

in = [Ufps 0 0 0; 

0 Ufps-Zalphadot 0 0; 

0 -l*Malphadot 1 0; 

0 0 0 1 ]; 

% Aircraft Matrix 

an = [Ufps*Xu Xalpha 0 -l*g*cos(thetanaut); 

Ufps*Zu Zalpha Ufps+Zq -l*g*sin(thetanaut); 
Ufps*Mu Malpha Mq 0; 

0 0 1 0 ]; 

% Control Matrix 

bn = [Xde Zde Mde 0]'; 

a = inv(in)*an; 
b = inv(in)*bn; 
c = eye(size(a)); 
d =[0 0 0 0]'; 

t = 0:01:100; 

[y,x,t] = step(a,b,c,d,l,t); 

xO = y(find(t— 15),:); 
xO = xO/xO(2)*5/57.3 
[y,x,t] = initial(a,b,c,d,x0,t); 

clg 

figure(l) 

%plotting alpha and theta, short period 
plot(t,y(:,2)) 

%title('Short Period Homogeneous Response'); 
grid 

axis([0 15 -.2 .2]); 
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gtext('alpha - rad'); 
hold on 

plot(t,y( : ,3));gtext('q - rad/sec'); 
xlabel('Time, sec'); ylabel('Response Gain'); 
hold off 
pause 

figure(2) 

%plotting u/U and theta, long period 
plot(t,y(:,l)) 

%title('Long Period Homogeneous Response'); 

grid 

axis; 

gtext('uAJ'); 
hold on 

plot(t,y(:,4));gtext('theta - rad'); 
xlabel('Time, sec'); ylabel('Response Gain'); 
hold off 
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5. doublet, m 



% This program uses EXPM function to find response to pitch doublet 

clear 

dderiv 

% Inertial Matrix 

in = [Ufps 0 0 0; 

0 Ufps-Zalphadot 0 0; 

0 -l*Malphadot 1 0; 

0 0 0 1 ]; 



% Aircraft Matrix 

an = [Ufps*Xu Xalpha 0 -l*g*cos(thetanaut); 
Ufps*Zu Zalpha Ufps+Zq -l*g*sin(thetanaut); 
Ufps*Mu Malpha Mq 0, 

0 0 1 0 ]; 

bn = [ Xde Zde Mde 0 ]'; 

a = inv(in)*an; 
b = inv(in)*bn; 



tl = 1.0; 


% start of doublet, sec 


t2 =2.0; 


% midpoint of doublet, sec 


t3 =3.0; 


% end of doublet, sec 


dO = -1; 


% unit elevau nput (1 rad) [t.e.u] 


tim = 15; 


% set end time 


t =0:0.05:tim; 



x = zeros(4,length(t)); % initialize solution matrices 
xl = zeros(4,length(t)); 
x2 = zeros(4,length(t)); 
x3 = zeros(4,length(t)); 

for i = l:length(t) 

ift(i) >=tl 
de(i) = dO; 
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xl(:,i) = dO*(expm(a*(t(i)-tl)) - eye(size(a)))*inv(a)*b, 
end 

ift(i) >= t2 
de(i) = de(i) - 2*d0; 

x2(:,i) = -2*d0*(expm(a*(t(i)-t2)) - eye(size(a)))*inv(a)*b; 
end 

ift(i) >= t3 
de(i) = de(i) + dO; 

x3(:,i) = dO*(expm(a*(t(i)-t3)) - eye(size(a)))*inv(a)*b; 
end 

x(:,i) = xl(:,i) + x2(:,i) + x3(:,i); 
end 

ymin = -2*abs(d0); 
ymax = 2*abs(d0); 

V = [0 tim ymin ymax]'; 

figure(l) 

plot(t,de);grid;axis(V) 

xlabel('Time, sec'); ylabel('Elevator, Rad'); 

%title('Elevator Input') 
pause; axis 

figure(2) 

%plotting alpha and pitch rate 
plot(t,x(2,:)) 

%title('Doublet Response') 

grid;gtext('alpha') 

hold on 

plot(t,x(3 , :));gtext('q') 

xlabel('Time, sec');ylabelCResponse Gain') 

hold off 
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6. spbode.m 

% This program uses the BODE function to find the 
% short-period response (transfer-function gain and phase) 
% to a harmonic input 

clear 

dderiv 

% Inertial Matrix 

in = [Ufps 0 0 0; 

0 Ufps-Zalphadot 0 0; 

0 -l*Malphadot 1 0; 

0 0 0 1 ]; 



% Aircraft Matrix 

an = [Ufps*Xu Xalpha 0 - 1 *g*cos(thetanaut); 

Ufps*Zu Zalpha Ufps+Zq -l*g*sin(thetanaut); 
Ufj)s*Mu Malpha Mq 0; 

0 0 1 0 ]; 

bn = [ Xde Zde Mde 0 ]’; 

a = inv(in)*an; 
b = inv(in)*bn; 
c = eye(size(a)); 
d = zeros(size(b)); 

w = logspace(- 1 , 2 , 1 50); 

[mag,phase,w] = bode(a,b,c,d,l,w); 

for i = 1 :length(w) 
for j = 1 :4 
if phase(i j) > 0 
phase(i j) = phase(i j) - 360; 
end 
end 
end 

clg 

figure(l) 

% Plotting alpha and pitch-rate gains 
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semilogx(w,mag( 1 ));grid;gtext(’alpha’) 
%title('Short-Period Frequency Response’) 
hold on 

semilogx(w,mag(:,2));gtext('q') 
xlabel('F requency, rad/sec');ylabel('Gain') 
hold off 
pause 

figure(2) 

% Plotting alpha and pitch-rate phase angles 
semilogx(w,phase( : , 1 ));grid;gtext('alpha') 
%title('Short-Period Frequency Response') 
hold on 

semilogx(w,phase(:,2)),gtext('q') 
xlabel('Frequency, rad/sec'),ylabel('Phase, deg') 
hold off 
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7. n bode.m 



% This program uses the BODE function to find 
% the normal-acceleration response 
% (transfer-function gain and phase) to a harmonic input 

clear 

dderiv 

% Inertial Matrix 

in = [Ufps 0 0 0; 

0 Ufps-Zalphadot 0 0; 

0 -l*Malphadot 1 0; 

0 0 0 1 ]; 

% Aircraft Matrix 

an = [Ufps* Xu Xalpha 0 -l*g*cos(thetanaut); 
Ufps*Zu Zalpha Ufps+Zq -l*g*sin(thetanaut); 
Ufps*Mu Malpha Mq 0; 

0 0 1 0 ]; 

% Control Matrix 

bn = [Xde Zde Mde 0]'; 

a = inv(in)*an; 
b = inv(in)*bn; 
c = [ Zu Zalpha Zq 0 ]; 
d = Zde; 

w = logspace(-l,2,150); 

[mag, phase, w] = bode(a,b,c,d,l,w); 

for i = 1 :length(w) 

mag(i) = mag(i)/g*pi/180; % converting to g's/deg 

if phase(i) > 0 
phase(i) = phase(i) - 360; 
end 
end 

% Plotting load factor/deg elevator input 
figure(l) 

semilogx( w, mag) ;grid 

xlabel(TYequency, rad/sec'); ylabel('Normal Acceleration Gain’) 
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%title('Normal Acceleration Response to Harmonic Input’) 
pause 

% Plotting normal acceleration phase angle 
figure(2) 

semilogx(w,phase);grid 

xlabel(’Frequency, rad/sec’); ylabel('Normal Acceleration Phase, deg’) 
%title(’Normal Acceleration Response to Harmonic Input 1 ) 
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B. LATERAL-DIRECTIONAL DYNAMICS 



1. lat_dir.m 

% Solves the full 4x4 lateral-directional response 

clear 

dderiv 

Idel latdir.dia 
diary latdir.dia 
in = [Ufps 0 0 0 

0 1 0 -1 *Ixz/Ixx 

0 0 10 

0 -l*Ixz/Izz 0 1]; 

an = [YB Yp g*cos(thetanaut) Yr-Ufps 



LB 


Lp 


0 


Lr 


0 


1 


0 


0 


NB 


Np 


0 


Nr]; 



a = inv(in)*an; 

[wn,zeta] = damp(a) 

[x,d] = eig(a) 
r= [0 0 0 0]; 

for i = 1 :4 
r (i) =d(i,i); 
end 

x2 = zeros(4,2); 
for j = 1:4 

x2(j,l) = abs(x(j,l)); 
x2(j,2) = 1 80/pi*angle(x(j, 1 )); 
end 

xx = x2(l,l); 
xy = x(3,2); 
xz = x(4,3); 
x3(:,l) = x2(:,l)/xx; 
x3(:,2) = x2(:,2) - x2(l,2); 
x3(:,3) = x(:,3)/xy; 
x3(:,4) = x(:,4)/xz; 
x3 

diary off 
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2. rudkick.m 



% Response to unit rudder impulse 

clear 

dderiv 

in = [Ufps 0 0 0 

0 1 0 -l*Ixz/Ixx 

0 0 10 

0 -l*Ixz/Izz 0 1]; 

% Plant Matrix 

an = [YB Yp g*cos(thetanaut) Yr-Ufps 



LB Lp 0 


Lr 


0 1 0 


0 


NB Np 0 


Nr]; 


% Control Matrix 
bn = [Ydr Yda 




Ldr Lda 
0 0 

Ndr Nda]; 




a = inv(in)*an; 
b = inv(in)*an; 
c = eye(size(a)); 





d = zeros(size(b)); 



t = 0:. 05:10; 

[x,y,t] = impulse(a,b,c,d, 1 ,t); 

% plotting sideslip (beta) and bank angle (phi) 
plot(t,y(:, 1 ));grid;gtext(’beta’) 
hold on 

plot(t,y( : , 3 ));gtext('phi') 

xlabel('Time, sec');ylabelCResponse Gain') 

%titleCRudder Kick Response') 

hold off 

pause 

% integrate p for final bank angle 
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phi = 0; 
t = 0:. 05:100; 

[x,y,t] = impulse(a,b,c,d,l,t); 

for i = l:length(t) 

phi = phi + y(i,2)*.05; 
end 

phi 
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3. sideslip, m 

% Solves the full 4x4 lateral-directional response for 
% initial sideslip condition 

clear 

dderiv 

in = [Ufps 0 0 0 

0 1 0 -l*Ixz/Ixx 

0 0 10 

0 -l*Ixz/Izz 0 1]; 

an = [YB Yp g*cos(thetanaut) Yr-Ufps 
LB Lp 0 Lr 

0 10 0 
NB Np 0 Nr]; 

bn = [ Ydr Yda 
Ldr Lda 
0 0 

Ndr Nda ]; 

a = inv(in)*an; 
b = inv(in)*bn; 
c = eye(size(a)); 
d = zeros(size(b)); 

t = 0:. 05:15; 



% calculate sideslip per bank angle 

f = [ 0 0 -1*CL ]'; 
k = [ CnB Cndr Cnda 
C1B Cldr Clda 
CyB Cydr Cyda]; 

perphi = inv(k)*f; 
betaperphi = perphi(l); 

PHI = 10*pi/180; 

BETA = betaperphi*PHI; 
xin = [ BETA 0 PHI 0 ]’ 
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[x,y,t] = initial(a,b,c,d,xin,t); 

plot(t,y(:, 1 ;grid ;gtext('bet a') 

%title('Dutch Roll Response for Sideslip I. C.s') 
hold on 

plot(t,y(:,3),'-');gtext('phi') 
hold off 

xlabel(’Time in Seconds') 
ylabel('Amplitude') 
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4. roll.m 



% Solves the full 4x4 lateral-directional for roll response 
% (transfer function gain and phase) due to 
% a harmonic aileron input 

clear 

dderiv 

in = [Ufps 0 0 0 

0 1 0 -1 *Ixz/Ixx 

0 0 10 

0 -l*Ixz/Izz 0 1] 

an = [YB Yp g*cos(thetanaut) Yr-Ufps 
LB Lp 0 Lr 

0 10 0 
NB Np 0 Nr] 

bn = [ Ydr Yda 
Ldr Lda 
0 0 

Ndr Nda ] 

a = inv(in)*an 
b = inv(in)*bn 
c = eye(size(a)) 
d = zeros(size(b)) 

w = logspace(-l,2,150); 

[mag,phase,w] = bode(a,b,c,d,2,w); 

for i = 1 :length(w) 
for j = 1:4 
if phase(ij) > 0 
phase(i,j) = phase(i j) - 360; 
end 
end 
end 

clg 

figure(l) 

% Plotting roll angle gain 
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semilogx(w, mag( : , 3 ));grid 
%title(’Harmonic Aileron Input') 
xlabel('Frequency, rad/sec');ylabel('Response Gain') 
pause 

figure(2) 

% Plotting roll angle phase 
semilogx(w, phase( : , 3 ));grid 
%title('Harmonic Aileron Input') 
xlabel('Frequency, rad/sec');ylabel('Phase in Degrees’) 
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